﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace CivLib.CivNum
{
  class NewtonRaphson : SolverBase
  {
    enum FiniteMethod { ForwardDifference, CentralDifference }
    public double m_currentSlope;
    public double m_currentError;
    public int m_varIndex = -1;
    double m_preturbation;
    FiniteMethod m_finiteDifferenceMethod = FiniteMethod.ForwardDifference;


    /// <summary>
    /// The preturbation used in the finite difference
    /// </summary>
    public double Preturbation
    {
      get { return m_preturbation; }
      set { m_preturbation = value > 0 ? value : m_preturbation; }
    }

    /// <summary>
    /// The type of difference method to use
    /// </summary>
    private FiniteMethod FiniteDifferenceMethod
    {
      get { return m_finiteDifferenceMethod; }
      set { m_finiteDifferenceMethod = value; }
    }

    public NewtonRaphson(ISoluableEqn a_variables, double a_preturbation = 0.01,
      double a_tolerance = 0.0001, double a_maxIter = 200)
      : base(a_variables, a_tolerance, a_maxIter)
    {
      m_preturbation = a_preturbation;
      m_tolerence = a_tolerance;
      m_maxIterations = a_maxIter;
    }

    void UpdateGradient()
    {
      double forward;
      switch (m_finiteDifferenceMethod)
      {
        case FiniteMethod.ForwardDifference:
          m_eqn[m_varIndex] += m_preturbation;
          forward = m_eqn.Residual();
          m_eqn[m_varIndex] -= m_preturbation;
          m_currentSlope = (forward - m_currentError) / m_preturbation;
          break;
        case FiniteMethod.CentralDifference:
          m_eqn[m_varIndex] += m_preturbation / 2;
          forward = m_eqn.Residual();
          m_eqn[m_varIndex] -= m_preturbation;
          m_currentSlope = (forward - m_eqn.Residual()) / m_preturbation;
          m_eqn[m_varIndex] += m_preturbation / 2;
          break;
      }

    }

    /// <summary>
    /// This method will go through the VarVec and find the unconstrained
    /// variable to solve for it.  If there is more than one unconstrained
    /// variable it will throw a SoverException.
    /// </summary>
    /// <returns>True if converged, false if not</returns>
    public override void Solve()
    {
      m_varIndex = -1;
      for (int i = 0; i < m_eqn.Count; i++)
      {
        if (!m_eqn.Constrained[i])
        {
          if (m_varIndex >= 0)
          {
            throw new SolverException(this.GetType().ToString() +
              ": both " + m_eqn.VarName[m_varIndex] + " and " +
              m_eqn.VarName[i] + " are unconstrained.");
          }
          m_varIndex = i;
        }
      }
      Solve(m_varIndex);
    }

    public void Solve(int a_varIndex)
    {
      m_varIndex = a_varIndex;
      for (int i = 0; i < MaxIterations; i++)
      {
        m_currentError = m_eqn.Residual();
        if (Math.Abs(m_currentError) < m_tolerence)
          return;
        UpdateGradient();
        double delta = m_currentError / m_currentSlope;
        m_eqn[m_varIndex] = m_eqn[m_varIndex] - delta;
      }
      throw new SolverException(this.GetType().ToString() +
        "did not converge, current error: " + m_currentError +
        " current slope: " + m_currentSlope);
    }

  }
}
